#include "../frame.h"
#include "../ini/velocity.h"
#include "../ini/position.h"
#include "../sys/rand.h"
#include "../io/xyz.h"
#include "../io/gro.h"
#include "../ini/crystal.h"
#include "../step/verlet.h"

#include <time.h>

#define ERROR program_result = -1; goto END

int main( int argc, char *argv[] )
{
	
	int program_result = 0;
	
	mds_frame    frame = MDS_FRAME_INIT;
	const size_t np    = argc > 1 ? strtol(argv[1], NULL, 10) : 100;
	const mds_vector box = { 10, 11, 12 };
	
	fprintf( stderr, "#particles= %u\n", (unsigned)np);
	mds_srand( time(NULL) );
	
	if( mds_frame_load( &frame, np, MDS_NUM_SETS_FOR_VERLET ) < 0 )
	{
		mds_perror( stderr, "mds_frame_load(%u,MDS_SETS_FOR_VERLET)\n", (unsigned)np);
		ERROR;
	}
	
	{
		mds_vector_set position = frame.sets[MDS_POSITION];
		mds_vector_set velocity = frame.sets[MDS_VELOCITY];
		register size_t i;
		
		for(  i = 0; i < frame.num_particles; ++i )
		{
			const Real mass = 12 + 10 * mds_alea();
			position[i].info.inv_mass = 1.0 /  mass;
			velocity[i].info.mass     = mass;
		}
		
		mds_velocity_initialize( position, velocity, np, 300.0 );
		mds_position_initialize( position, np, mds_pbc_xyz, &box);
		
		mds_xyz_save("rbox.xyz", position, np, "random box", 0);
		mds_gro_save("rbox.gro", "random box", position, velocity, np, &box, 0);
		
	}
	
	{
		FILE *fp = fopen("anint.txt", "wt" );
		if( fp )
		{
			Real x;
			for( x = -4.1; x <= 4.1; x += 0.05 )
			{
				fprintf( fp, "%g %g\n", x, mds_anint(x) );
			}
			fclose( fp );
		}
		fprintf( stderr, "anint(-0.49) = %g\n", mds_anint( -0.49 ) );
		fprintf( stderr, "anint(-0.51) = %g\n", mds_anint( -0.51 ) );
		
	}
	
	
	
END:
	MDS_MEMORY();
	mds_frame_quit( &frame );
	MDS_MEMORY();
	
	return program_result;
	
}
